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Abstract 

P 

The  objective  of  this  research to  compare  a  quasi- 
analytical,  potential  flow/three-degree-of-freedom  model  to 
an  implicit-Euler  algorithm  for  the  calculation  of  store 
trajectories.  The  implicit  algorithm  uses  a  cell-centered/ 
finite-volume,  spatial  discretization  applied  to  the  Euler 
equations,  written  in  time-dependent,  curvilinear- 
coordinates.  A  flux-differencing  Roe  scheme  is  employed  to 
find  the  split-fluxes  and  the  Steger/Warming  flux-vector 
method  is  used  to  calculate  the  f lux-Jacobians .  The 
potential  flow,  and  implicit-Euler  algorithm  are  combined 
with  a  three-degree-of-freedom  algorithm  to  evaluate  the 
planar,  f reef all  trajectories  of  a  simple  store  shape.  The 
research  uses  two  different  grid-modification  techniques  in 
the  implicit  algorithm  evaluation. 

Data  collected  for  both  grids^utilizea  the  minimum 
time-step  in  the  three-degree-of-freedom  algorithm  for  a 
Courant  number  of  10.  Two  test  cases  involved  updating  the 
flux-Jacobians  after  every  time-step  and  only  once  during 
every  1000  iterations.  The  effect  of  multiplying  the 
minimum  time-step  by  factors  of  2,  4,  6,  8,  10,  and  100  were 
also  examined. 

The  potential  flow  and  implicit  algorithm  trajectories 
didn't  compare  very  closely.  The  various  At  and  Jacobian- 
update  results  matched  rather  closely. 


COMPUTATION  OF  PLANAR  STORE  TRAJECTORIES 
USING  AN  ADAPTIVE  GRID  PROCEDURE 


Chapter  I.  Introduction 


Weapon  separation  testing  is  an  integral  part  of  the 
test  and  evaluation  carried  out  on  every  United  States  Air 
Force  Tactical  and  Strategic  aircraft.  Prior  to  flight 
test,  the  only  reliable  method  of  determining  the  trajectory 
of  a  weapon  is  v;ind-tunnel  testing  assisted  by  computer 
analysis  (25:4.34).  Accurate  and  robust  Computational  Fluid 
Dynamics  (CFD)  algorithms  would  make  additional  weapon 
trajectory  data  readily  available. 

Wind-tunnel  testing  is  a  very  costly,  time  consuming, 
and  rigid  process.  V?ind-tunnel  models  take  several  weeks  to 
fabricate  and  are  very  expensive.  If  modifications  to  the 
aircraft  or  v/eapon  occur  during  the  testing  process,  wind- 
tunnel  data  may  be  invalidated.  Wind-tunnel  tests  generally 
take  several  weeks  to  plan  and  prepare,  and  once  completed 
the  data  are  only  valid  for  the  weapon  and  aircraft 
configurations  tested.  Flight  test  planning  uses  these  data 
in  determining  the  best  and  safest,  separation-flight  test. 

Flight  testing  is  also  costly;  it  is  hazardous  as  well. 
Before  flight  testing  occurs,  both  the  engineer  and  pilot 
need  an  idea  of  the  trajectory  of  the  weapon.  Only  wind- 
tunnel  data  provides  this  trajectory  information.  If 
modifications  to  the  weapon  or  aircraft  occur,  the  v;ind- 
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tunnel  test  data  may  be  invalid,  sometimes  requiring  further 
wind-tunnel  testing  or  a  more  drawn  out  and  hazardous 
flight-test  program.  If,  however,  the  weapon  and  aircraft 
could  be  computationally  modeled,  changes  to  the  aircraft  or 
weapon  could  be  easily  incorporated  and  new  trajectory  data 
computed  relatively  quickly. 

There  is  still  much  research  required  before  CFD 
algorithms  will  be  able  to  provide  accurate  weapon 
separation  data.  This  research  is  an  attempt  in  furthering 
the  development  of  a  CFD  algorithm  to  the  weapon  separation 
problem.  I  will  investigate  the  problem  of  computationally 
moving  a  simple  geometric  shape  (ellipse)  through  a  two- 
dimensional  flow-field  using  Whitfield's,  dynamic-grid, 

Euler  algorithm  (28) .  The  investigation  will  be  limited  to 
planar  trajectories  and  will  ignore  viscous  effects.  The 
work  presented  will  exercise  the  Whitfield  algorithm  and 
compare  its  results  against  a  quasi-analytical,  potential- 
flow  result. 

The  organization  of  this  thesis  is  as  follows.  Chapter 

II  presents  the  potential-flow  solution  and  differential 
form  of  Euler  equations.  The  transformation  of  the  Euler 
equations  from  rectangular  Cartesian-coordinates  to  general 
time-dependent,  curvilinear-coordinates  is  explained.  This 
transformation  allows  dynamic  grids  to  be  utilized.  Chapter 

III  presents  the  implicit,  upwind,  finite-volume,  flux- 
vector  scheme.  This  chapter  also  discusses  the  use  of  flux- 
differencing  with  Roe-averaging  and  the  use  of  flux-limiters 
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to  achieve  schemes  with  higher-order  accuracy.  Chapter  IV 
gives  a  brief  discussion  of  the  equations  of  motion  used  to 
generate  the  Three-Degree-of-Freedom  (3-DOF)  algorithm.  The 
grid-manipulation  technique  used  in  this  research  is  also 
presented.  Chapter  V  presents  the  results  of  the  research. 
The  research  compares  the  potential-flow  results  with 
results  from  the  Whitfield  algorithm.  It  continues  by 
examining  the  effect  of  modifying  the  flux-Jacobian  update 
methodology  and  the  minimum  time-step  on  the  Whitfield 
algorithm.  Chapter  VI  closes  with  conclusions  and 
recommendations  on  further  research. 

I . 1  Summary  of  Current  Knowledge 

Wind-tunnel  testing  with  computer  simulation  is 
currently  the  only  method  of  resolving  the  trajectory  of  a 
weapon  from  an  aircraft  without  flight  testing.  Wind-tunnel 
testing  obtains  the  forces  and  moments  on  the  weapon  used  in 
a  Six-Degree-of-Freedom  (6-DOF)  program  utilized  by  the 
3246th  Test  Wing/TYES  to  compute  store  trajectories  (17:1). 
The  6-DOF  program  determines  the  trajectory  of  a  weapon  very 
accurately  when  compared  with  flight-test  data.  This 
technique  has  been  reliably  used  for  many  years.  The 
trajectories,  however,  are  limited  to  the  aircraft  and 
weapon  configurations  and  flight-test  conditions  performed 
during  wind-tunnel  testing.  Only  changes  in  the  moments-of- 
inertia,  center-of-gravity  (COG) ,  and  weight  of  the  weapon 
can  be  accounted  for  in  the  6-DOF  algorithm.  Physical 
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modifications  to  the  weapon  or  aircraft  require  new 
aerodynamic  data  to  account  for  changes  in  the  flow-field 
about  the  weapon  and  aircraft.  These  data  are  in  turn  used 
to  calculate  the  correct  forces  and  moments  present  on  the 
weapon  and  its  trajectory. 

The  wind-tunnel  test  limitations  could  be  overcome  by 
the  development  of  CFD  algorithms,  which  could  accurately 
model  weapon  separation.  The  Wright  Laboratory,  Armament 
Directorate,  Eglin  Air  Force  Base  is  developing  CFD 
algorithms  to  model  the  separation  of  weapons  from  aircraft. 
The  Whitfield  algorithm  was  developed  by  Mississippi  State 
University  under  contract  as  part  of  this  work  for  the 
Armament  Directorate  (6:1).  The  algorithm  accurately  models 
the  forces  and  moments  on  the  weapon  generated  by  the 
surrounding  flow-field.  In  References  1,  2,  6  through  8, 
and  27  through  29,  algorithm  results  have  shown  good 
agreement  with  wind-tunnel  data  for  various  store  and 
airfoil  shapes.  These  results  contained  data  for  both 
stable  and  dynamic  grids  in  both  planar  and  three- 
dimensional  cases.  The  algorithm  can  account  for  viscous 
and  compressibility  effects  in  subsonic,  transonic,  and 
supersonic  flow.  This  algorithm  has  the  added  ability  of 
utilizing  blocked  grids  (8:31).  This  allows  a  tailored  grid 
to  be  generated  around  each  object,  such  as  an  aircraft  and 
weapon,  which  is  to  be  modelled.  These  tailored  grids  can 
be  coupled  and  a  complete  solution  of  the  flow-field 
calculated.  For  example,  the  grid  around  the  weapon 
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requires  modification  as  the  weapon  trajectory  is  determined 
by  the  aerodynamic  forces  and  launch  conditions. 

Figure  1  shows  an  example  of  how  the  grid  can  be 
modified  as  a  store  (an  ellipse  for  this  example)  is 
dropped.  The  Whitfield  algorithm  has,  however,  only  been 
examined  for  dynamic  grids  with  specified  trajectories 
(7:13).  These  specified  trajectories  move  the  entire  grid. 
The  algorithm  hasn't  been  tested  for  modified  grids  similar 
to  that  shown  in  Figure  1 


(a)  (b) 

Figure  1:  Grid  at  t  =  1  (a)  and  Grid  at  t  =  2  (b) 

My  research  extends  the  testing  of  the  Whitfield 
algorithm  into  the  realm  of  unspecified  trajectories.  This 
requires  modification  of  the  grid  after  each  time-step  to 
account  for  the  new  position  of  the  ellipse  as  specified  by 
the  3-DOF  model.  The  calculated  trajectory  will  be  compared 
with  trajectories  calculated  using  the  quasi-analytical. 
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potential-flow  solution.  This  assumes  a  good  comparison 
exists  between  low-speed,  compressible  flow  and  quasi- 
analytical,  potential  flow.  By  adding  the  3-DOF  algorithm, 
and  eventually  the  6-DOF  algorithm,  to  the  Whitfield 
algorithm,  a  true  CFD  algorithm  could  result  that  would  fill 
a  void  in  the  current  capabilities  of  weapon  test  and 
evaluation. 
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Chapter  II.  Theory 


II. 1  Potential-Flow  Solution 

The  initial  model  developed  for  this  research  uses  an 
ellipse  to  represent  a  weapon  and  an  analytical,  potential- 
flow  solution.  The  comparison  with  the  Whitfield  algorithm 
solutions  utilizes  the  analytical  solution  for  potential 
flow  about  an  ellipse  in  determining  a  stable  trajectory.  A 
stable  trajectory  is  defined  as  one  where  the  pitching 
motion  of  the  store  decreases  as  the  store  transitions  into 
f reef all.  Stable  stores  generally  require  the  COG  to  be 
forward  of  the  center-of-pressure.  The  potential-flow 
assumption  initially  limits  the  research  to  cases  of 
inviscid  and  incompressible  flow.  These  assumptions  are 
used  for  verifying  the  trajectory  algorithm. 

Lamb  presents  a  potential-flow  solution  for  the  flow- 
field  about  a  rotating  ellipse.  He  used  the  complex 
potential  that  describes  potential  flow  past  a  circular 
cylinder.  Lamb  transformed  this  flow,  using  complex 
variables,  into  flow  about  an  ellipse  (14:89).  Lamb 
presents  the  following  form  for  the  stream  function 

T  =  ^  exp(-5)  (Ub  sinr\  -  Va  cost]) 

+  -ici)(a+b)^  exp(-2^)  cos2r] 

+  c(U  sinhJi  sim\  -  V  coshJi  cost])  (1) 
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In  Eq.  1,  U  and  V  represent  the  two  freest ream,  velocity- 
components;  0)  is  the  angular  velocity  of  the  ellipse, 
measured  about  its  center;  a  and  b  are  the  semi-major  and 
semi-minor  axis  of  the  ellipse.  Figure  2  shows  the 
orientation  of  U,  V,  o),  a,  and  b. 


Figure  2:  Velocity  Component  Relationship 


c  is  a  parameter  defined  by 

c  H 


(2) 


The  semi-minor  axis  for  the  ellipse  used  in  this  research  is 
0.2a. 

^  and  T|  are  the  elliptic  coordinates  used  in  the 
complex  variable  transformation.  The  transformation 
equations  are 


X  =  C  cosh^  COSTl 

(3a) 

y  =  c  sinh^  sinx\ 

(3b) 

The  calculations  of  lift,  drag,  and  moment  on  the 
ellipse  use  a  trapezoidal  integration  of  the  surface 
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pressure.  The  calculation  assumes  no  atmospheric  pressure 
gradient.  The  surface-velocity  components,  u  and  v,  are 
calculated  from  the  stream  function,  Eq.  1,  by 


{4a) 


V 


dx 


(4b) 


or,  in  transformed  coordinates 


_  1 
dx  J 


dy 


i5a) 

(5b) 


The.  Jacobian,  J,  in  Eqs.  5  is 

J  =  x^y^-x^y^  (6) 

In  this  potential-flow  solution,  there  is  no 
circulation  about  the  ellipse.  Thus,  the  resulting  lift 
component  is  zero.  The  drag  is  also  zero,  since  the 
pressure  distribution  is  exactly  balanced.  The  calculated 
results  helped  validate  the  trapezoidal-integration 
technique.  The  moment  calculation  is  also  utilized  as  a 
comparison  of  the  integration  accuracy,  Batchelor  gives  an 
analytical  solution  for  the  moment  about  a  non-rotating 
ellipse  at  a  specified  angle  of  attack  (4:435), 
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M  =  -p%uv(a^  -  h^) 


(7) 


This  equation  was  used  to  verify  the  moment  calculation. 
Appendix  A  contains  results  of  the  trapezoidal-integration 
analysis . 


II. 2  The  Euler  Equations 

The  Euler  equations  are  a  simplified  frem  of  the 
Navier-Stokes  equations  describing  the  behavior  c:.  continuum 
fluid  flows.  The  Euler  equations  neglect  the  effects  of 
viscosity,  heat  transfer,  and  body  forces.  The  prime 
denotes  dimensional  variables  in  the  following  equations. 

The  conservative,  differential  form  of  the  Euler  equations 
in  Cartesian-coordinates  are  (29:1) 


+  df  ^  dg'  dh>  _  ^ 

dt'  dx'  dy'  dz' 


(8) 


v/here  q' ,  £' ,  g' ,  and  h'  are  defined  by 


qr'  =  Jp',  pv,  pv,  e'] 

T 

(9a) 

~  [p'u’ ,  p'u^  +  P',  p'v^u' ,  p^/u^. 

{e'+P')  u']  ^ 

(9b) 

-  [p'v' ,  p'u'v'.  p'v^ ,  +  P' ,  p'v/v^ , 

{e'+P')  v']  ’’ 

(9c) 

h'  =  [p'i/,  p'u^w' ,  p*v'v£ ,  p'vP'  +  P' , 

(e'+pO  t/jr 

(9d) 
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The  Cartesian  (x',  y'/  z')  velocity-components  are 
represented  by  {u' ,  v' ,  w') .  P'  and  p'  represent  the 
pressure  and  density,  respectively.  The  total  energy  of  the 
flow  e' ,  is  defined  through  the  perfect  gas  law  assumption 
as 

e'  =  +  — p'(u^  +  +  w^)  (10) 

y-l  2 


where,  y  is  the  ratio  of  specific  heats  (y  =  1.4).  More 
often,  the  Euler  equa*"  ions  are  written  in  nondimensional 
form.  The  nondimensional  variables  used  in  this  analysis 
are  defined  as  (2:7) 


(11a) 


(11b) 


(11c) 


where 


is  the  freestream  speed  of  sound.  The  variable,  1', 
represents  any  convenient  reference  length  (2:8)  .  The 
reference  length  is  2a  for  the  ellipse, 

Eq.  8  can  be  rewritten  in  nondimensional  form  as 


dt  dx  dy  dz 


(13) 


II. 2.1  Generalized  Curvilinear-Coordinate 
Transformation 

The  development  of  the  Euler  equations  used  the 
Cartesian-coordinate  system  (x,  y,  z)  .  In  most 
applications,  however,  this  coordinate  system  doesn't 
facilitate  the  accurate  numerical  treatment  of  boundary 
conditions.  It  is  also  insufficient  in  representing  complex 
geometries,  such  as  airfoils,  aircraft,  weapons,  etc. 
Instead,  a  boundary-conforming,  coordinate  system  should  be 
incorporated.  The  curvilinear-coordinates  are  explicit 
functions  of  the  Cartesian-coordinates 


5=5  {x,y,  z,  t) 

(14a) 

n  =  n  ix,y,  z,  t) 

(14b) 

C  =  C  (x,y,  z,  t) 

(14c) 

X  =  t 

(14d) 

These  curvilinear-coordinates  change  the  nondimensional  form 
of  the  Euler  equations  to 
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(15) 


dx  dK  dr\  dC~ 


where 


p 

pu 

pu 

puu  +  5^ 

pv 

F  =  J 

PVU  +  lyP 

pw 

pwu  + 

e. 

(e+P)  u  -  5c-P 

(16a) 


pr 

pw 

PUV  +  TljpP 

puW  + 

pvV  +  r\yP 

II 

pvf«r  + 

pwV  +  Ti^P 

pwjv  + 

(e+P)  V  -  x\^P 

(e+P)  W  -  CtP 

(16b) 


U,  V,  and  Wf  are  the  contravariant  velocities/  in  the  x\, 
and  ^-coordinate  direction.  They  are  defined  as 

(7  =  5c  ^•y'^ 

fvT  =  C(.  +  C_^u  +  CyV  +  (17c) 


The  Jacobian,  J,  represents  the  volume  of  each  three- 
dimensional  cell  in  physical  space  and  is 


The  metric  terms  are 
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=  J'^iz^y^  -  y^z^) 

(19a) 

Tly  =  tj"^  (x^z^  -  Z^Xf) 

(19b) 

=  J'''  (-^nn  " 

tlz  =  -  x^y^} 

(19c) 

5c  =  (-x,5^-y,5y-2x5P 

Tic  =  (-x,ii^-y,tiy-z,tiP 

(19d) 

Cx  = 

. (19e) 

(iy  =  -  z^x^) 

(19f) 

Cz  = 

(19g) 

Cc  =  (--^tCx  -  ^tCy  -  ^tCz) 

(19h) 

Janus  details  these  transformations  (13:92-94), 
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Chapter  III.  Algorithm  Development 


III. 1  The  Implicit  Algorithm 

The  algorithm  that  solves  the  Euler  equations  is  a 
combination  of  the  Roe  scheme  and  the  Steger/V7arming  flux- 
vector  splitting  scheme  (2:35) .  These  schemes  take 
advantage  of  characteristic  theory  as  applied  to  the 
hyperbolic  Euler  equations.  The  advantages  to  this  method 
are  the  elimination  of  numerical  viscosity  or  smoothing 
terms  and  an  increase  in  the  convergence  rate.  This 
algorithm  e:<ploits  the  accuracy  of  the  Roe  scheme  in 
defining  discontinuities  and  the  relative  computational  ease 
of  flux-vector  splitting  techniques  to  give  accurate 
solutions.  The  order  of  accuracy  of  the  algorithm  is  user 
specified  by  the  limiter  applied.  Limiters  and  their  order 
of  accuracy  are  discussed  in  Chapter  III. 3. 

III. 1.1  Flux-Vector  Splitting 

Flux-vector  splitting  developed  by  Steger  and  Warming 
(Reference  24)  in  the  1980s  has  been  widely  used  in  the 
numerical  solution  of  the  Euler  equations.  The  procedure 
accounts  for  the  direction  of  information  travel  determined 
from  characteristic  theory  to  compute  fluxes.  The  three- 
dimensional  Euler  equations  (Eq.  15)  are  a  system  of  five 
equations  that  have  five  characteristic  velocities  in  each 
of  the  three  coordinate  directions.  These  characteristic 
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velocities  are  resolved  from  the  quasilinear  form  of  Eq.  15 
(29:2) , 


^  4-  B  ^  -  0 


(20) 


where  the  flux- Jacobian  matrices,  A,  B,  and  C,  are 


A  =  ^ 

do 


B  = 


do 


C  = 


dH 

do 


(21) 


The  eigenvalues  of  A  are  the  characteristic  velocities  in 
the  ^-direction.  Similarly,  the  eigenvalues  of  B  and  C  are 
the  characteristic  velocities  in  the  T|-  and  ^-directions, 
respectively. 

F,  G,  and  H  are  identical  except  ^  appears  in  F,  T]  in 
G,  and  ^  in  H.  We  can  simplify  the  analysis  by  letting  K 
represent  either  F,  G,  or  H  (28:3).  We  can  then  define 


do 


(22) 


K  corresponds  to  A,  B,  and  C  depending  on  the  meaning  of  K. 
The  eigenvalues  of  the  matrix  K  are  (13:18) 


+  kyV  +  k^w  +  iCj.  =  P;t 


4  =  Pic  +  =  pjt  +  +  kl  +  kl)  2 


(23a) 

(23b) 


=  p;c  -  cm 


(23c) 
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where  here  c  is  the  speed  of  sound  and  k  is  either  T|,  or 
^  corresponding  to  A,  B,  or  C,  respectively. 

The  flux-vector,  K,  can  be  split  into  three  parts,  each 
part  corresponding  to  the  distinct  eigenvalues  of  K  given  in 
Eqs.  23.  Janus  gives  the  details  of  this  splitting  (13:98- 
100) .  The  flux-vector,  K,  is  then  written 

K  =  klK^  +  +  kUs  (24) 


where 


P 

P 

pu 

pu+pck^ 

pv 

pv+pcky 

Y 

pw 

*  2y 

pw+pck^ 

£  (u^  +  V^+K^) 
2 

e+p+pcd^ 

K.  = 


J 

2y 


P 

pu-pck^ 

pv-pcky 

pw-pck^ 

e+p~pcQ^ 


(25) 


and 


(26) 


=  k^u  +  kyV  +  k^w 


(27) 
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The  sign  of  X  in  Eqs.  23  controls  the  direction  that 
information  is  used  in  determining  the  analogous  portion  of 
the  flux,  Kj^,  v/here  i  =  1,  4,  and  5  in  Eqs.  25.  This  flux- 
vector  split  is  not  unique  and  other  techniques  for 
splitting  into  positive  and  negative  parts  are  possible 
(2:15).  The  splitting  used  by  Steger  and  Warming  is  given 
by 


II 

+  F" 

(28a) 

G  =  G* 

+  G- 

(28b) 

'H  =  H* 

+  H- 

(28c) 

Here,  the  plus  superscript  signifies  the  portion  of  the 
flux-vector  associated  with  the  non-negative  eigenvalues. 
The  minus  superscript  signifies  the  portion  of  the  flux- 
vector  associated  with  the  non-positive  eigenvalues 
(24:269) . 

III. 1.2  Finite-Volume  Formulation 

This  algorithm  uses  a  finite-volume  formulation,  thus 
allowing  for  the  modeling  of  arbitrary  configurations.  It 
can  better  manage  the  complicated  grid-structure,  than  the 
finite-difference  formulation  (2:15) .  The  computational 
domain  is  divided  up  into  many  small  "non-overlapping" 
cells.  In  this  way  Eq.  8,  can  be  evaluated  over  each  cell 
separately  or  over  the  entire  computational  domain.  The 
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resulting  solution  will  fulfill  the  Euler  equations  over  a 
specified  volume. 

The  dependent  variables  in  these  equations  are  stored 
at 'each  grid  point  or  cell  center.  The  flux-vectors  used  in 
this  finite-volume  formulation,  however,  are  required  at 
cell  surfaces.  Figure  3  shows  a  typical  portion  of  the 
computational  domain  and  its  nomenclature. 


V 


(i.  j+1/2) 

u 

O  ( 

(i.j) 

)(i+l/Z.  i) 

Figure  3:  Indexing  of  Cell  Centers  and  Surfaces 


One  discretized  form  of  the  three-dimensional, 
unsteady,  Euler  equations  (Eq.  15)  is 

At  A5  Ati 


+ 


=  0 


(29) 
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Eq.  29  can  be  written  more  compactly  as 


AO"  +  =  0  (30) 


where 


Ac?"  =  C?"*^  -  C?" 


(31) 


and 


»  pn*l  -  ''^*1/2, j.k  ~  ^-1/2. j.k) 


(32) 


The  central  differences  of  G  and  H  are  defined  similarly. 

The  n  superscript  indicates  values  at  the  current  time-step, 
while  the  n+1  superscript  indicates  unknown  values  at  the 
next  time-step. 

The  implicit,  flux-split  discretization  of  the  Euler 
equations  is  formed  by  substituting  Eqs.  28  into  Eq.  30  to 
give, 

AC>"  +  AT[fi{(F"+F-)"^^  +  a^(G*+G-)"*^  +  5j(ff"+H-)«*i]  =0  (33) 

The  flux-vectors  and  are  nonlinear 

functions  of  the  dependent  variables,  A  linearization 

procedure  using  a  local  Taylor  series  expansion  of  the  flux- 
vectors  about  the  time-step,  n,  was  suggested  by  Beam  and 
Warming  (5:89).  The  Taylor  series  expansion  is 
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-  0") 

+  0(At2) 

(34a) 

(F-)n*i  =  F-‘^+  A- 

-  0”) 

+  0(At2) 

(34b) 

where 


A* 


(35) 


Similar  linearizations  are  used  for  the  split  versions  of 

and  Substitution  of  Eqs.  34  into  Eq.  33  yields, 

[I  +  AtCSjA*-  +  SjA"-  +  +  8fC-*)]AC>"  =  -Ari?" 

(36a) 

where 

E"  =  [65  (F*  +  F-)n  +  +  G")"  +  +  H')”]  (36b) 

k”  is  called  the  residual.  The  dot  in  the  above  equation 
symbolizes  that  the  difference  operators  apply  to  the 
product  of  the  flux- Jacobian  matrices,  A,  B,  and  C,  with 
Ao”.  It  is  important  to  note  that  the  linearization  in  Eqs. 
34  is  second-order  accurate  in  time,  so  that  the  algorithm, 
as  represented  by  Eqs.  36  remains  second-order  accurate. 

The  f lux-Jacobians  on  the  left-hand-side  of  Eq.  36a  are 
determined  using  the  flux-vector,  split  scheme  developed  by 
Steger  and  Warming  (24) .  The  residual  term  implements  the 
Roe  Scheme  in  determining  the  values  of  the  split  fluxes. 
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The  f lux-Jacobians  would  be  computationally  very  costly  to 
obtain  using  the  Roe  Scheme.  Whitfield  states,  "in  all 
results  obtained  thus  far,  the  Jacobian  matrices 
corresponding  to  the  flux-vector  split  scheme  gave  improved 
convergence  rates  and  a  more  robust  scheme  than  did  the  Roe 
matrices.  The  reason  for  this  is  unclear  (28:21)."  This 
combination  of  the  Steger/Warming  and  Roe  schemes  were 
demonstrated  in  Reference  2  to  provide  accurate,  dynamic- 
grid  results  compared  with  wind-tunnel  data. 

III. 1.3  Factorization 

A  numerical  solution  of  Eqs.  36  at  this  point  cannot  be 
•realistically  obtained  because  of  the  large-banded  matrix  of 
the  system.  The  left-hand-side  of  the  equation  requires  the 
inversion  of  a  large  matrix  and  an  exact  solution  would  be 
too  numerically  intensive.  A  way  to  alleviate  this  problem 
is  to  factor  the  operator  on  the  left-hand-side  of  Eq.  36a 
into  several  operators.  Whitfield  in  Reference  27  and 
Anderson  in  Reference  1  have  applied  many  such 
factorizations  to  Eqs.  36.  These  schemes  include  a  six, 
three,  and  two-factor  scheme.  The  two-factor  scheme  used  in 
this  research,  separates  the  implicit  operator  into  one 
factor  containing  all  the  flux-Jacobians  with  the  non- 
negative  eigenvalues  and  the  other  factor  containing  all  the 
Jacobians  with  the  non-positive  eigenvalues.  The  two-factor 
scheme  can  be  written  as 
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+  +  [J+At(8^A‘*  +  5^S*-  +  5jC--)  ]  Af?"=-ATi2' 


This  equation  can  be  solved  in  the  following  steps 


[J  +  AtCS^A"*  +  b^B*'  +  5(C"-)]A0*  =  -AtJ?' 


(38a) 


[I  +  A'c(85A--  +  6^B--  +  5^C--)]ACJ"  =  AO* 


(38b) 


=  0"  +  Af?” 


(38c) 


The  spatial-differences  on  the  left-hand-side  of  Eq.  37  are 
only  required  to  be  first-order  accurate  to  retain  second- 
order  accuracy  overall.  The  implicit  algorithm  is  trying  to 
drive  A(3  0  in  the  steady-state  calculations,  so  higher- 

order,  spatial-difference  terms  are  not  necessary  on  a  fine 
grid.  The  spatial-difference  terms  are  evaluated  using  one- 
point  extrapolation  as  given  by  (2:21) 


(39a) 


b^(A  AQ")^  -  [A'  AQi^j^k  ~ 


(39b) 


The  metrics  at  appropriate  cell  face' s  are  required  to 
evaluate  the  Jacobian  matrices.  The  A‘*'(Q^,-  >,)  and 

terms  use  the  metrics  at  the  (i+l/2,j, /c)  and 
{i-l/2,j,k)  cell  faces,  respectively.  The  spatial- 
differences,  5^(B'^AQ),  5^(B“AQ),  5^(C'‘’AC?),  and  5^(C“A(})  are 
similarly  defined. 


The  two-factor  scheme  requires  only  the  solution  of  a 
sparse,  block  lower  and  upper-triangular  matrix  at  each 
time-step  (2:21).  The  solution  to  Eq.  38a  is  calculated  by 
a  simple-forward  substitution  and  the  solution  of  Eq,  38b  by 
a  simple-backward  substitution.  This  approach,  however, 
requires  the  storage  of  three,  large,  f lux-Jacobians 
matrices  at  each  of  the  cell  centers  at  any  given  time- 
level.  Anderson  presents  the  stability  analysis  for  the  two 
factor  scheme  in  Reference  1,  His  analysis  suggests  Courant 
number  doesn't  affect  the  two-factor  scheme  and  it  retains 
stability  for  Courant  numbers  up  to  35  (1:4). 

III. 2  Flux-Difference  Splitting 

Thus  far,  the  discussion  has  focused  on  the  Steger  and 
Warming  flux-vector,  splitting  scheme,  which  is  used  to 
evaluate  the  left-hand-side  of  Eqs.  36.  Whitfield 
(Reference  28)  also  incorporates  the  flux-difference  scheme 
developed  by  Roe,  in  the  solution  to  the  right-hand-side  of 
Eqs.  36.  The  flux-difference  scheme  is  incorporated,  since 
it  has  been  shown  to  capture  shocks  and  contact 
discontinuities  in  as  few  as  one  cell.  Whitfield's 
algorithm  was  designed  to  investigate  a  wide  range  of  fluid 
flow  problems,  which  has  been  enhanced  by  the  incorporation 
of  the  Roe  scheme.  Bram  van  Leer  demonstrated  a  first- 
order,  Roe  scheme  required  only  one-fourth  to  one-half  as 
many  grid  points  to  accurately  capture  shocks  and  contact 
discontinuities  as  a  second-order,  flux-vector  split  scheme 


24 


(2:23)  .  These  advantages  are  being  exploited  by  the 
incorporation  of  the  flux-difference  scheme  into  the  flux- 
vector  scheme. 


III. 2.1  Riemann  Problem 

The  very  essence  of  the  flux-difference  splitting 
scheme  is  the  solution  of  the  local  Riemann  problem  (28:3). 
Considering  the  initial-value  problem  for  a  hyperbolic 
conservative  system  of  equations  written 


dt 


(40) 


g(x,  0)  =  g(x)  for  -<»  <  x  <  <» 


(41) 


where  q{x,t)  is  a  column  vector  of  m  unknowns  and  f(u),  the 
flux,  is  a  vector-valued  function  of  m  unknowns.  The 
finite-volume  discretization  of  Eq.  40  is  given  as 

Qi '  -  Qi  ^  _  Q  (42) 

At  Ax 


The  term  represents  the  average  value  at  time,  =  nAt, 
in  the  interval  (i-l/2)Ax  <  x  <  (i+l/2)Ax.  The  numerical 
flux-vector,  i+\/2i  evaluated  at  the  grid-cell 
interface. 

Godunov  developed  a  procedure  to  advance  the  solution, 
qj_,  to  the  next  time-level  by  solving  a  set  of  Riemann 
problems  at  each  cell  interface  (2:24) .  The  Riemann  problem 
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can  be  strictly  applied  to  a  set  of  conservative  equations 
if  initial  data  are  given  over  semi-infinite  intervals.  For 
example/  consider  the  initial  data  (Eq.  41)  given  for  the 
initial-value  problem  (Eq.  40) 


<7(x 


,0)  =  [; 


Ql 

<Jr 


for  x<  (i+l/2)Ax| 
for  x>  (i+l/2)Axj 


(43) 


where  qj^  =  and  q^  =  q^+i  for  all  cell  interfaces, 

X  =  (i+l/2)Ax.  Godunov  utilized  an  integral  average  over 
the  interval,  (i-l/2)Ax  <  x  <  (i+l/2)Ax  to  solve  the  Riemann 
problem  at  the  cell  interfaces  to  obtain  (2:25).  The 

time-interval  selected  in  Godunov' s  scheme  must  ensure  no 
interaction  between  neighboring  Riemann  problems  occurs, 
resulting  in  an  exact  solution  to  the  Riemann  problem  over 
the  complete  domain.  Approximate  Riemann  solvers  have  been 
developed  to  reduce  the  complex  structure,  but  retain  the 
important  properties  of  the  solution. 

III. 2. 2  Roe's  Approximate  Riemann  Solver 

Roe  developed  a  method  utilizing  the  Riemann  solver  for 
solving  a  set  of  linear-conservation  laws  (2:26) .  He 
computed  the  exact  solution  of  the  Riemann  problem  for  the 
linear-hyperbolic  system,  instead  of  finding  an  approximate 
solution  to  the  complete  equations.  The  linearized  system 
is  a  rewritten  form  of  Eqs.  40-41, 
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(44) 


dt: 


+  =  0 


Qix,0) 


for  X  <  o\ 
q^  for  X  >  Qj 


(45) 


where  Aiqj^,qp)  is  a  constant  matrix.  The  constant  matrix  is 
chosen  to  represent  local  interface  conditions.  Roe 
required  the  matrix,  Aiq^^^q^) ,  to  have  the  following 
properties  (2:26): 

(i) 

(ii) 

(iii) 

(iv) 

Once  the  A  matrix  is  assembled,  its  eigenvalues  can  be 
estimated  as  the  wavespeeds  of  the  Riemann  problem. 
Eigenvalues  of  the  matrix  A  are  represented  by  .  The 
eigenvalues  are  distinct  and  arranged  in  increasing  order 
(i.e.  <  ...  <  X^) .  The  corresponding  left  and  right 

eigenvectors  are  represented  by  and  r^ ,  respectively.  It 
can  be  shown  that  the  change  in  the  dependent  variable  is 
proportional  to  the  right  eigenvectors  of  the  flux-Jacobian 
matrix,  whether  the  matrix  ^  is  constant  or  not  (2:27) . 

This  means  dq  across  each  characteristic  curve,  connected 


It  forms  a  linear  mapping  from  the  vector  space, 
q,  to  the  vector  space,  f 

As  ->  g,  A(qj^rqp)  approaches  A(g),  where 

A  =  Bf/dq 

For  any  g^,  q^:  A(qj^,q^)  *  (g^^  -  q^)  =  f{q^)  -  f{q^) 
The  eigenvalues  of  A  are  linearly  independent. 
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with  a  certain  eigenvalue/  is  proportional  to  the  right 
eigenvector,  associated  with  that  eigenvalue.  Thus, 

Qr-  Ql= 

An  intermediate  state,  q^,,  between  any  two  characteristic 
curves,  K  and  K+1,  can  be  computed  in  terms  of  the  right 
eigenvectors  of  A  in  the  following  way 

The  interface  differential,  dq,  is  proportional  to  the  right 
eigenvector  of  the  A  matrix  and  considering  that  the  right 
eigenvectors  are  calculated  from  (A  -  X^I)r^  =  0,  the 
following  equation  can  be  developed  (2:27) 

Adq^^^  -  (48) 

Therefore,  the  total  change  across  all  characteristic  curves 
is  calculated  by 

-  Qi)  =  ^a^X^r^  =  fjj  -  f i  (49) 


The  quantity,  dfj  =  ajX^r^,  represents  the  change  in  the 
flux-vector  across  the  characteristic  curve-  The  flux- 
vector  is  required  at  each  cell  interface  to  solve  the 
hyperbolic  system  of  conservation  laws.  Let's  say,  a  cell 
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interface  is  located  at  x  =  0;  f'i+1/2  is  the  numerical 
flux-vector  at  that  interface/  which  can  be  computed  by 
either  of  the  following  two  expressions 


(50a) 


fL-,„  =  fu-  taA*^r^ 


(50b) 


where  -  and  +  represent  the  summation  over  the  negative  and 
positive  eigenvalues  of  “wavespeeds"/  respectively  (2:28). 
The  average  of  the  two  flux-vector  representation  is 


fi.1/2  =  I  4  *  4  -  paj\XJ\zi 


where  is  the  magnitude  of  the  jth  wave.  Whitfield,  in 
Reference  28,  has  shown  that  by  the  hyperbolic  nature  of  the 
system,  the  matrix  A  can  be  diagonalized  through  a 
similarity  transformation  as  A  =  TAT~^ .  The  diagonal 
matrix.  A,  is  composed  of  the  eigenvalues  of  A.  T  is  a 
matrix  whose  columns  are  the  right  eigenvalues  of  A  and  T~^ 
is  the  matrix  whose  rows  are  the  left  eigenvalues  of  A.  Eq. 
51  can  be  written  as 


*  fjr  -  |A(gt,<3rR)  -  gj] 


where 
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j'A  CTjj)  I  ~  ~  ^  Qr) 


(53) 


and 


=  T 


m 


y-i 


(54) 


An  averaging  process  of  the  dependent  variables  is 
heeded  in  order  to  malce  the  A  matrix  satisfy  Roe's 
requirements  as  cited  in  Chapter  II I. 2. 2.  The  Roe-averaged 
variables  for  the  three-dimensional  Euler  equations 
incorporate  information  from  either  side  of  the  cell 
interface  as  denoted  by  the  subscripts  L  and  R,  (28,10) 


P  = 


(55) 


u  = 


Pc  ***  Pr 


pI'^ 


Pr 


(56) 


V  = 


^l/2, 


,1/2, 


Pl  Vj  ^  PR  ^R 

^1/2  .  .1/2 
Pl  +  Pj? 


(57) 


w  = 


(58) 


H  = 


PR^f^R 


pr 


(59) 


a2  =  =  (Y-1) 

P 


jy  - 


u^+v^+w^ 


(60) 
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H  =  —  (e+P) 

P 


(61) 


e  =  — ^  +  -I  (u2+v2+ivr2)  (g2) 

y-l  2 

where  P  is  the  pressure;  p  is  the  density;  u,  v,  and  w  are 
the  velocity  components;  e  is  the  total  energy;  and  H  is  the 
total  enthalpy. 

The  one-dimensional  analysis  can  be  extended  to  three- 
dimensions  and  used  with  the  transformed  curvilinear- 
coordinate  system  following  the  appropriate  transformations 
(Eqs.  14) .  The  wave  can  be  split  into  the  three 
curvilinear-coordinate  directions  and  the  associated 
eigenvalues  of  the  Euler  equations  determine  the  wave 
motions  normal  to  the  cell  interfaces.  The  approximate 
Riemann  solver  is  then  applied  at  each  cell  interface 
(2:30). 

III. 3  Higher-Order  Spatial  Accuracy 

The  numerical  solution  in  Eq.  52  yields  a  first-order 
accurate  scheme,  which  we  would  like  to  improve  to  second- 
or  third-order  accuracy.  We  should  simultaneously,  however, 
maintain  the  desirable  sharp,  shock-capturing 
characteristics  without  introducing  any  pre-shock 
oscillations.  A  higher-order  accurate,  total-variation- 
diminishing  (TVD)  scheme  using  Roe-averaging  was  introduced 
by  Osher  and  Chakravarthy  and  has  been  adapted  for  the 
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present  work  by  Arabshahi  (2:31)  .  A  thorough  discussion  of 
the  implementation  of  these  higher-order  methods  is  in 
Reference  24.  A  brief  discussion,  however  will  be  provided 
here. 

A  third-order  accurate,  numerical  flux-vector  at  the 
cell  interface,  i+1/2,  can  be  obtained  by  the  addition  of 
correction  terms  to  the  first-order  accurate  flux  in  Eq.  52. 
This  calculation  uses  Roe-averaged  variables  and  metric 
terms  for  the  computation  of  eigenvalues  and  eigenvectors  at 
cell  interfaces.  The  numerical  flux-vector  is  given  as 

+  t  (-1,1)  -  l;  (3 , 1) 

+  S1^[l;(-1,-1)  -  Lj(l,3)]r^/2  (63) 

where 

and 

~  ^+1/2  *  "  Oi-x) 

®J,i*l/2  “  ^i+l/2  *  ((?i+i  “  Oi) 

1*3/2  ~  ^1*1/2  *  (C*l+2  “  ^1*0 


(65a) 

(65b) 

(65c) 
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When  computing  higher-order-accurate/  upwind-biased  schemes, 
overshoot  and  undershoot  are  expected  in  the  vicinity  of  the 
shock.  A  limiter,  jt  such  as  the  one  in  Eg.  63,  can  be 
implemented  to  reduce  the  scheme  to  a  fully  one-sided  scheme 
for  first-order  accuracy  in  the  shock  region  (2:31).  This 
helps  to  eliminate  the  overshoot  trait  of  higher-order 
schemes . 

Within  the  Whitfield  algorithm,  three  limiting 
operations  are  available.  They  are  the  minimum-modulus 
(minmod) ,  superbee,  and  van  Leer  limiter.  Each  of  these  was 
used  in  the  research. 

The  minmod  limiter  uses  Eg.  63,  where  the  limiting 
operator,  is  defined  by 

Lj(l,n)  =  minmod  (oJ.i+1/2,  (66) 

The  operator  "minmod”  is  given  as 

minmod[x,y]  =  signi^x)  *  max(0,  min[|x|,  y  sign{x)]  )  (67) 

The  parameter,  b,  is  called  the  compression  parameter  and  is 
calculated  from  the  function  (28:16) 

b  =  (68) 

1-ijr 


where  y  is  the  accuracy  parameter.  The  minmod  limiter 
contains  two  arguments,  x  and  y.  When  the  two  arguments  x, 
and  y,  are  of  opposite  sign,  the  operator  value  is  zero. 
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When  the  arguments  are  of  the  same  sign,  the  operator 
chooses  the  argument  with  the  smaller  absolute  value  (2:33). 

The  subscript,  i+1/2,  in  Eq.  63  indicates  where  the 
Roe-averaged  variables  and  metrics  are  calculated.  The 
eigenvalues,  and  left  and  right  eigenvectors,  and 

r-^,  are  evaluated  using  Roe-averaged  variables.  The  flux- 
vectors,  f(Qj)  and  in  Eq.  63  are  computed  using  the 

dependent  variables.  The  principle  part  of  the  truncation 
error  of  this  higher-order  scheme  can  be  shown  to  be  (2:32) 


The  accuracy  of  the  scheme  can  be  altered  merely  by  varying 
the  accuracy  parameter,  y.  If  \j/  =  1/3,  a  TVD  scheme  of 
third-order  accuracy  results.  A  choice  of  \|f  =  -1  gives  a 
second-order  upwind  TVD  scheme  (28:16). 

The  superbee  limiter  is  another  limiter  used  in  this 
algorithm.  The  limiting  operator,  jt  this  scheme  is 

Lj  (C,  n]  =  cmplim  [o5,i+,/2,  <Jj,i*n/2] 


where 


cinplim\.x,y]  =  sign[x]  * 


max  (  0,  min[|x|,py  sign  (x)],  min(p|x(,y  sign{x)]  )  (71) 


34 


Here,  (3  is  the  compression  parameter.  It  is  generally 
between  1  <  P  <  2,  and  normally  P  =  2  is  used  (28:17) . 

The  van  Leer  limiter  is  the  third  limiter  used  in  the 
Whitfield  algorithm.  Here,  the  limiter  is 

=  vanlim 


where 

vanllm  U,y]  =  ^  *  l-^^l  (73) 

x+y 

The  superbee  and  van  Leer  limiters  are  independent  of  the 
accuracy  parameter,  \}/. 

The  flux-difference  scheme,  with  the  limiters  just 
discussed,  is  used  to  evaluate  the  residual  on  the  right- 
hand-side  of  Eq.  37  in  the  Whitfield  algorithm.  The  left- 
hand-side  of  Eq.  37  is  calculated  using  the  first-order, 
upwind,  flux-vector  splitting  scheme.  The  accuracy  of  the 
solution  is  very  good  and  the  convergence  rate  for  steady- 
state  improves  by  combining  the  two  schemes  (2:35).  When 
the  algorithm  is  evaluating  a  dynamic  problem  the  higher- 
order  accuracy  terms  remain  applicable.  Further  explanation 
of  this  algorithm  can  be  found  in  Reference  2. 

III. 4  Boundary  Conditions 

The  boundary  conditions  are  very  critical  in  the 
determination  of  the  ultimate  form  of  the  flow-field. 
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Whitfield  concluded  the  use  of  characteristic-variable 
boundary  conditions  produced  the  best  results  (2:36),  The 
Euler  equations  must  be  put  into  characteristic  form  before 
determining  the  boundary  conditions.  The  derivation  of  the 
Characteristic  form  of  the  Euler  equations  can  be  found  in 
Reference  29.  There  are  three  boundary  conditions  required 
for  this  investigation:  subsonic  inflow,  outflow,  and 
impermeable  surface.  Only  the  surface  boundary  condition  is 
dependent  bn  whether  the  grid  is  stationary  or  dynamic. 

III. 4.1  Subsonic  Inflow 

The  subsonic  inflow  case  is  characterized  by  the  first 
four  eigenvalues  being  positive  and  the  fifth  is  negative 
with  flow  in  the  increasing  computational  ^c-coordinate 
directions,  T),  or  The  first  three  and  fifth 
eigenvalues  are  negative  and  the  fourth  is  positive  for  flow 
in  the  decreasing  computational  /c-coordinate  direction 
(29:7).  The  subsonic  inflow  conditions  are  (29:8) 

■Pi  =  ^ ± Po<^o (“a" “i)  (74a) 

Pi  =  Pa  +  (74b) 

cl 

Uj,=  u^±  (74^,) 

*  ^  ^  PoCo 
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(Jc^  kl)  2 


where  i  is  x,  y,  or  z,  depending  on  which  computational 
coordinate  is  being  evaluated.  The  point  a  is  outside  the 
computational  domain,  point  Jb  is  on  the  computational  domain 
and  point  I  is  inside  the  computational  domain  (29:8). 

Figure  4  shows  the  relationship  between  these  points. 


IXOir  FLOW 


Flcrjr  in  +X  Direction  How  in  --K  Direction 

Figure  4:  Computational-Coordinate  Schematic  for 
Boundary  Conditions 


The  plus/minus  sign  option  corresponds  to  the  sign  of  the 


first  three  eigenvalues. 


III. 4. 2  Subsonic  Outflow 

The  subsonic  outflow  conditions  are 


Pb  =  Pi 


(76a) 


p,  =  p, . 


(76b) 


“Jb  =  “a  ± 


Pa  -  Pb 

Po^o 


(76c) 


V.  =  ±  £  l^P. 

^  PoCo 


(76d) 


Wu  =  tvr 


2j  -  vva  = 


Po^o 


(76e) 


where  the  point  a  is  inside  the  computational  domain  and 
point  I  is  outside  the  computational  domain  (29:9).  The 
plus/minus  signs  have  the  same  meaning  as  in  Eqs.  74. 

I I I . 4 . 3  Impermeable  Surface 
The  impermeable  surface  boundary  condition  is 
characterized  by  the  first  three  eigenvalues  being  zero,  the 
fourth  is  positive  and  the  fifth  is  negative  (13:43).  The 
boundary  conditions  for  a  stationary,  impermeable  surface 
are 

Pb=  Px  ^  Po^o  *  ^z^x'> 
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P.  =  p.  *  7' 

(77b) 

cl 

Uy  -  +  KyVj.  + 

(77c) 

Vj.  -  Ky 

(77d) 

Wj.  -  Kx  *  ^z^z^ 

{lie) 

where  the  subscript  r  refers  to  the  center  of  the  first  cell 
inside  the  boundary  or  the  reference  value.  In  Eq.  77a,  the 
minus  sign  is  used  if  r  is  in  the  positive  >c-direction  from 
the  boundary.  The  plus  sign  is  used  if  r  is  in  the  negative 
ic-direction  from  the  boundary  (29:9). 

The  dynamic  impermeable  surface  condition  varies  for 
the  stationary  condition  due  to  the  addition  of  the  time 
term.  Eqs.  77  become  (7:12) 

-  Pb  ~  Pr  ~  ^  (78a) 


“6 


=  Pr  - 

Ap 

cl 

(78b) 

Uz  -  Kx 

Ap 

Po^o 

(78c) 

Vz  -  Ry 

Ap 

Po^o 

(78d) 

'^z  - 

Ap 

n  n 

(78e) 
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When  the  dynamic  grid  is  implemented,  the  metrics  and  flux- 
Jacobians  are  normally  updated  after  every  time-step  to 
account  for  the  motion  of  the  store. 

The  last  concept  to  be  addressed  is  phantom  points. 
Phantom  points  can  lie  inside  the  boundary  of  the  ellipse  or 
outside  the  computational  domain.  Figure  5  shows  the 
relationship  between  the  boundaries,  phantom  points,  and 
reference  points. 


\  I 


Boundary  Surface 


Figure  5:  Reference  and  Phantom  Point  Schematic 


Points  a  and  1  use  the  phantom  point  concept  in  the 
determination  their  values  (13:45).  The  subscript  p 
represents  the  phantom  point  and  the  values  for  the  phantom 
points  are  calculated  using 
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(79) 


Op  =  20i,  -  0^ 


where  0  can  be  p,  p,  u,  v,  or  w.  The  i-subscript  refers  to 
the  first  cell  inside  the  computational  domain  and  can  be  a, 
1,  or  r  depending  on  the  boundary  condition.  For  example, 
the  phantom  points  for  the  impermeable  surface  would  be 
calculated  using  (29:9) 


Ij.  T  2PoCo  +  iCyVj.  +  K^Wj.) 

(80a) 

rp  r  r  2 

(80b) 

Co 

“r  -  2iC;,  (ic^U^  +  £yV^  + 

(80c) 

-  2Ky 

(80d) 

Wr  -  2ic^  {Rx^r  +  RyVj.  +  RxWj) 

(80e) 

III. 5  Planar  Application  of  the  Three-Dimensional 
Whitfield  Algorithm 

The  three-dimensional  Whitfield  algorithm  has  been 
adapted  for  application  to  this  planar  problem.  This  was 
accomplished  by  generating  a  one-cell,  three-dimensional 
grid,  where  the  x-  and  y-locations  of  the  grid  points  in 
each  z-plane  were  identical.  The  two  x-y  planes  where 
spaced  a  unit  distance  apart  in  the  z-direction.  The  z- 
component  of  velocity,  w,  was  also  set  to  zero.  This 
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allowed  the  ellipse  to  be  treated  as  a  two-dimensional  body 
in  the  3-DOF  algorithm. 


42 


Chapter  IV.  Trajectory  Algorithm  and  Grid  Modification 


IV. 1  Equations  of  Motion 

A  3-DOF  model  was  added  to  the  potential  flow  and 
Whitfield  algorithms.  The  equations  of  motion  for  the 
ellipse  are  derived  from  Newton's  second-lav/-of-motion, 
which  states; 

A  body  acted  upon  by  a  force  moves  in  such  a  manner 
that  the  time  rate  of  change  of  momentum  equals  the 
force  (11:90). 

The  rates  of  change  are  with  respect  to  an  inertial  space. 
This  law  can  be  written  in  two  vector  equations; 

Zf  =  (81) 

dt 

^  (82) 

dt 

where  m  is  the  mass  of  the  body,  0  is  the  velocity  vector, 
and  P  is  the  angular-momentum  vector.  Assuming  that  the 
earth  is  flat  and  that  gyroscopic  effects  are  negligible, 
Eqs.  81  and  82  reduce  to  a  system  of  three  equations  of 
motion  (10:32): 

-m  ^ 

=L-W^  cosy  (84) 
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(85) 


=  M 


where  y^QQg  is  the  speed  of  the  COG;  D  is  the  dra  •  I.  is  the 
lift;  M  is  the  moment  about  the  COG;  Wj.  is  the  weight  of  the 
ellipse;  is  the  moment-of-inertia  in  the  z-direction;  y 
is  the  angle  between  inertial  x-axis  and  the  velocity 
vector,  V’^og'  ®  angle  between  the  inertia  axis 

and  the  body  axis.  The  lift  and  drag  used  is  calculated  for 
the  wind  axis  system  which  is  parallel  to  the  inertial  axis 
system.  Figure  6  shows  the  relationships  between  the 
different  angles  and  the  orientation  of  the  ellipse  between 
the  inertial  and  body  axes. 


Yb 


Figure  6:  Orientation  of  Ellipse  with  3-DOF  Angles 
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9  and  y  are  measured  positive  counter-clockwise  and  a  is 
measured  positive  clockwise. 

Eqs.  83-85  represent  an  initial-value  problem  involving 
two  first-order  and  one  second-order  nonlinear,  ordinary- 
differential-equations.  This  problem  can  be  reduced  to 
three  first-order,  ordinary-differential-equations  as  shown 
in  Gallaway  by  (10:33); 

"  (8«> 


where. 


dt 


Q 


(87) 


By  using  Eqs.  86  and  87,  Eqs.  83-85  are  altered  to  give. 


dg  _  M 
dt  ■  I, 

dVcog  ^  ^  siny) 

dt  m 

dy  _  (L  -  cosy) 
~dt  " 


(88) 


(89) 


(90) 


These  three  differential  equations  can  be  integrated 
using  a  variety  of  techniques.  In  Callaway's  dissertation, 
both  Euler's  explicit  method  and  Adam-Bashforth  formulae  are 
used  to  solve  Eqs.  88-90  (10:35).  Pauletti's  6-DOF 
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algorithm  incorporates  a  Runge-Kutta  integration  scheme  to 
solve  the  equations  of  motion  (20:22).  A  6-DOF  algorithm  is 
utilized  by  the  3246th  TW/TY  to  determine  store  trajectories 
using  wind-tunnel  data.  This  algorithm  exercises  either  the 
Runge-Kutta  or  the  Adam-Moulton  integration  formulae,  which 
can  be  selected  by  the  user  (17:31) . 

Euler's  explicit  method  was  selected  for  use  in  this 
research  for  its  simplicity.  Euler's  method  is  only  first- 
order  accurate.  Callaway  found  no  significant  difference  in 
results  using  the  Euler  or  Adam-Bashforth  methods  in  his 
research  (10:35).  Euler's  method  leaves  us  with  four 
equations  (10:34) 


= 

^  cog 


V‘ 


cog 


at 


(91) 


yU+l  =  y"  +  At-^" 

at 


(92) 


=  gfl  +  A  t-f?” 

ot 


(93) 


0n*l  =  0n  +  At 


6^^ 

dt 


(94) 


Results  from  Eqs.  88-90  and  Eqs.  86  and  87  are 
used  to  calculate  and  Eq.  94  uses 

the  results  from  Eq.  93  to  calculate  0^’*’^,  utilizing  the 
relationship  in  Eq.  87.  Also,  the  follov/ing  relationship  is 
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used  in  the  potential-flow  algorithm  to  update  the  value  of 
(X  at  the  new  time-step  (10:29)  . 


a 


tan"^ 


V-,+  V  IcosyI 
^cofflsinvl 


(95) 


The  relative  velocity  seen  by  the  ellipse  is  the  difference 
between  the  ellipse  orientation  angle,  0,  and  the  velocity 
angle  of  the  COG,  y.  The  new  parameters  calculated  in 
Eqs.  91-95  are  used  in  the  next  iteration  to  calculate  the 
forces  and  moment  on  the  ellipse. 

The  3-DOF  algorithm  was  developed,  so  the  trajectory  of 
the  store  is  measured  relative  to  an  aircraft  flying 
straight  and  level  at  constant  velocity.  The  algorithm  was 
developed  in  this  way,  so  it  could  later  be  incorporated 
into  a  complete  aircraft/store  configuration.  The  initial 
conditions  are  constant  velocity  with  no  angular  velocity. 
The  store  is  then  allowed  to  freefall. 

For  the  quasi-analytical,  potential-flow  algorithm,  the 
lift,  drag,  and  moment  values  are  calculated.  They  are  then 
used  in  Eqs.  88-90  to  calculate,  dq^/dt,  dy’^/dt, 
and  89^/9t.  dq^/dt  =  (ip  in  Eq.  1.  These  results  are  used 
in  Eqs.  91-93  to  calculate  and  Using 

Eq.  87  and  is  calculated  from  Eq.  94.  Eq.  95  is 

then  used  to  calculate  the  value  of  and  are 

used  to  compute  the  new  location  and  orientation  of  the 
ellipse  relative  to  the  constant  velocity  aircraft.  The 
process  is  repeated  with  the  new  values  of  and 
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Qjn+l  used  in  Eq.  1  to  calculate  the  values  of 

lift,  drag,  and  moment  at  the  next  time-step. 

The  3-DOF  algorithm  is  similarly  integrated  into  the 

implicit  Euler  algorithm.  The  lift,  drag,  and  moment  values 

are  calculated  by  the  implicit  algorithm  and  and 

0n+l 

are  calculated  using  the  same  above  procedure, 

Y^^^,  and  are  then  used  to  re-orient  the  ellipse  and 

the  grid.  The  grid-manipulation  procedure  is  discussed  in 
depth  in  a  later  section. 

IV. 2  Center-of-Gravitv  and  Moment-of-Inertia 

The  COG  and  moment-of-inertia  calculations  for  the 
ellipse  are  required  to  create  a  stable  trajectory.  In 
general,  a  stable  trajectory  requires  the  COG  to  be  forward 
of  the  center-of-pressure  (26) .  The  ellipse  was  assumed  to 
be  composed  of  two  different  materials  having  different 
densities.  The  materials  are  linearly  distributed  as  shown 
in  Figure  7 . 


y 


Figure  7:  Mass  Distribution 
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Pj_  and  p2  are  constants.  This  mass  distribution  keeps  the 
GOG  in  the  y-axis  constant  at  zero. 


The  GOG  was  calculated  using 


where 

dn?!  =  2p^y(x)  dx 
dm2  ~  2p2y(x)  dx 


(96) 


(97) 

(98) 


y(x)  is  defined  by  the  equation  for  an  ellipse 


(99) 


where  a  and  b  are  the  semi-major  and  semi-minor  axes, 
respectively.  Therefore,  y(x)  is 

-  (100) 

a2  } 


y{x)  =  \b^ 


/ 

=  i 
\ 


Eq.  96  can  be  integrated  to  determine  the  GOG  location, 
depending  on  the  value  of  X2  selected. 

The  moment-of-inertia  is  calculated  using  (11:204) : 
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(101) 


+  y(x).^)din,  +  +  y{x)^)din2 


The  ellipse  used  in  this  research  was  composed  of  two 
materials  with  p-j^  =  2700.0  kg/m^,  the  density  of  aluminum, 
and  =  11342.0  kg/m^,  the  density  of  iron. 

IV. 3  Grid  Manipulation 

In  this  investigation,  the  grid  is  updated  to  reflect 
the  ellipse's  new  location  as  it  moves  through  an  arbitrary 
trajectory.  Grid  manipulation  was  not  required  for  the 
quasi-analytical,  potential-flow  algorithm,  as  explained  in 
Chapter  IV. 2. 

Two  different  grid-manipulation  techniques  were  used. 
The  first  techniques  moves  the  entire  grid  through  the 
arbitrary  trajectory  calculated  by  the  3-DOF  algorithm. 

This  is  the  same  procedure  used  in  References  2,  6,  and  8. 
Figure  8  shows  an  exaggerated  example  of  how  the  entire  grid 
translates  with  the-  ellipse.  The  first  technique  will  be 
referred  for  the  remainder  of  the  text  as  Grid  1 . 

The  second  technique  utilizes  an  elliptical  grid- 
generator  in  determining  the  new  grid  after  each  time-step. 
The  grid-generation  algorithm  was  developed  by  AFWAL/FIMM 
(David  J.  Amdahl)  and  has  been  modified  to  run  on  the 
STELLAR  ST-2000  by  Capt  Mark  Driver.  The  program  uses 
Thompson's  techniques  for  generating  grids  describe  in 
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Figure  8:  Example  of  Grid  1  Trajectory 


The  ellipse  is  moved  the  distance  calculated  within  the 
grid.  The  outer  boundary  of  the  grid  remains  stationary  and 
the  branch-cut  in  the  0-grid  is  moved  using  a  geometric 
progression  routine.  The  branch-cut  in  the  0-grid  is  where 


Figure  9:  Conversion  of  Computational  Grid  to  Physical  Grid 


Once,  the  new  location  of  the  ellipse  and  branch-cut 
are  calculated  the  elliptic  grid-generation  algorithm 
regenerates  the  grid.  The  ellipse  starts  in  the  position 
similar  to  that  in  Figure  9c.  Figure  10a  shows  an  example 
of  the  ellipse  and  branch-cut  moved  as  a  result  of  the  3-DOF 
algorithm.  Figure  10b  shows  the  grid  after  it  has  been 
modified  by  the  elliptic  grid-generation  algorithm.  This 
grid-manipulation  technique  will  be  referred  to  as  Grid  2 
for  the  remainder  of  the  text.  The  discontinuity  in  the 
cells  on  either  side  of  the  branch-cut  is  a  result  of  the 
elliptic  grid-generation  algorithm.  The  discontinuity  could 
not  be  removed  v;ith  the  grid-generation  algorithm. 


(a)  (b) 

Figure  10:  Generation  Procedure  for  Grid  2 

By  using  the  elliptic  grid-generator,  the  uniform  outer 
boundary  can  remain  stationary.  This  allows  for  the 
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utilization  of  a  block  grid-structure  in  the  CFD  algorithm 
to  support  an  aircraft/store  separation  problems.  The 
Whitfield  algorithm  was  developed  for  complex  blocked  grid- 
structures  (8:32).  The  algorithm  allows  for  the  generation 
of  a  grid  around  a  specific  aircraft  with  "holes"  in  the 
grid  where  stores  could  be  added  to  complete  the  physical 
domain.  The  only  requirement  is  the  location  of  the  grid 
points  along  the  block  boundary  must  be  coincident.  A  grid- 
generation  routine  similar  to  the  Grid  2  technique  could  be 
used  to  keep  the  block-grid  boundary  stationary;  yet  still 
move  the  store  through  an  arbitrary  trajectory.  Belk 
describes  the  block-grid  technique  in  his  dissertation 
(8:32-53) . 
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Chapter  V:  Results 


This  section  presents  the  data  gathered  for  the  quasi- 
analytical/  potential-flov/  algorithm  and  two  grid- 
modification  cases.  The  preliminary  investigation  discusses 
the  initial  analyses  used  in  the  determination  of  the  Mach 
number  and  limiter  used  in  the  research.  From  this 
investigation,  no  limiter  was  utilized  in  the  actual 
trajectory  calculations.  Mach  =  0.3  was  selected  for  use  in 
the  implicit  Euler  algorithm.  A  comparison  between  the 
trajectory  trends  of  the  potential  flow  and  implicit 
algorithm,  for  Grid  1  and  2,  was  performed.  In  all  cases, 
the  ellipse  pitched  nose-down  and  then  begin  to  oscillate 
with  varying  periods.  The  variations  depended  on  time-step 
used  and  the  algorithm  and  grid  type.  The  robustness  of  the 
algorithm  and  solution  accuracy  were  also  investigated  in 
light  of  time-step  variation  and  changes  in  the  Jacobian- 
update  methodology.  Time-step  variation  influenced  the 
pitching  motion  of  the  ellipse  more  than  the  trajectory  of 
the  COG.  The  Jacobian-update  methodology,  however,  did  not 
effect  the  COG  trajectory  or  pitching  motion  of  the  ellipse. 

V.l  Preliminary  Analysis 

The  analysis  began  with  a  comparison  of  the  coefficient 
of  drag,  C^,  values  versus  grid  refinement  for  a  stationary 
ellipse.  An  exact  solution  to  the  Euler  equations  would 
result  in  no  Cj^  or  Cq[  calculated  about  a  symmetric  body. 
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Any  variation  is  a  result  of  either  numerical  viscosity  or 
numerical  error  introduced  by  the  algorithm.  In  the  initial 
results,  negative  occurred  using  the  limiters  in  the 

Whitfield  algorithm.  was  used  in  determining  the  grid 
refinement  necessary  in  generating  accurate  results.  The 
initial,  grid-refinement  investigation  used  the  Van  Leer 
limiter  in  calculating  the  C^.  The  maximum  grid-refinement 
was  a  201x51  grid  with  a  \|r  =  -1.  The  began  to  decrease 
through  zero  as  shown  in  Figure  11. 


Figure  11:  Cd  vs  Grid  Refinement  (with  Limiter) 

The  minmod  and  superbee  limiters  were  also  investigated 
attempting  to  resolve  the  negative  problem.  The  minmod 
limiter  generated  a  second  order  (MINMODl,  y  =  -1)  and  third 
order  (MINMOD2,  y  =  1/3)  solution.  The  superbee  limiter 
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generated  a  second-order  solution  with  a  \j/  =  0.  Figure  11 
also  shows  these  values. 

Coefficient  of  pressure,  Cp,  curves  are  shown  for  the 
55x51  grid.  Figure  12  shows  the  variations  in  C„  with  each 
limiter.  The  potential  flow  and  no  limiter  case  are  shown 
for  comparison. 

The  specific  reason  for  these  negative  drag  values 
cannot  be  explained  by  this  author.  A  numerical  error  is 
added  to  the  implicit  algorithm  by  these  limiters.  The 
numerical  error  negatively  impacts  the  Cp  on  the  ellipse, 
which  is  used  in  the  calculation  of  coefficient  of  lift,  Cj^ 
and  C^.  These  results  are  not  reasonable  compared  to  what 
nature  tells  us  about  drag.  Therefore,  limiters  were  not 
used  in  the  trajectory  investigation. 

Pulliam  hints  at  the  negative  drag  results  in  his  AIAA 
article  (Reference  21) ,  but  he  also  has  no  specific 
explanation  for  its  occurrence  (22) .  An  explanation  for  the 
negative  drag  coefficient  and  the  non-zero  lift  solutions  as 
discussed  by  Pulliam  requires  further  research.  Multiple 
telephone  conversations  with  the  developers  of  the  implicit 
algorithm  did  not  bring  to  light  any  additional  information 
into  the  reason  for  the  negative  drag  values  (3,23). 

By  using  no  limiter  in  the  implicit  algorithm,  the 
negative  C^  problem  is  alleviated.  Figure  13  shows  the  C^j^ 
versus  grid-refinement  plot.  A  201x201x2  grid  was  initially 
chosen  for  the  research.  These  data  was  calculated  using 
the  Cray  2  supercomputer. 
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Figure  12:  Cp  vs  X  (For  Different  Limiter) 
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By  evaluating  the  coefficient  of  pressure,  Cp,  plots 
for  the  ellipse,  the  Mach  number  used  in  the  research  was 
determined.  From  these  data,  the  Mach  number  chosen  for  the 
research  was  0.3.  This  Mach  number  was  selected  because  of 
its  relatively  smooth  Cp  curve.  Mach  =  0.3  is  also  within 
the  Mach  range  where  flow  can  still  be  considered 
incompressible.  Therefore,  the  implicit  algorithm  results 
should  more  closely  match  those  of  quasi-analytical, 
potential-flow  results.  Appendix  B  details  these  results. 

A  problem  with  computer  availability  required  the 
remainder  of  the  research  to  be  performed  on  the  Stellar  ST- 
2000  computer.  Memory  requirements  inherent  to  the  implicit 
algorithm  required  the  grid  refinement  to  be  reduced  to 
101x101x2 . 
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V.1.1  Preliminary  Grid  Analysis 


Two  different  grids  were  used  during  the  evaluation  of 
the  implicit  algorithm.  Both  grids  have  an  outer  radius  of 
10  chords.  Grid  1  is  a  simple  0-grid  generated  using 
Eqs.  3.  Grid  2  was  generated  by  taking  Grid  1  and 
manipulating  it  using  the  elliptic  grid-generator  discussed 
in  Chapter  IV. 3.  The  steady-state  solution  for  Grid  1  and 
Grid  2  was  calculated  using  the  implicit  algorithm.  The  p- 
residual  {r’^)  converged  to  a  value  of  1.375x10“®  for  Grid  1 
and  2.579x10“®  for  Grid  2,  the  asymptotic,  steady-state 
value.  This  was  the  starting  point  for  the  trajectories 
gathered  in  this  research.  Figures  14  and  15  show  a  close- 
up  of  Grid  1  and  Grid  2  near  the  ellipse. 


2.5 


Y  0.0 


-2.5 

-2.5  0.0  2.5 

X 

Figure  14:  Grid  1,  101  x  101 
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X 

Figure  15:  Grid  2,  101  x  101 

The  implicit  solutions  from  both  grids  gave  slightly 
different  values  for  and  The  ellipse  in  Grid  1  had  a 

=  0.0000  and  a  =  0.1006,  whereas  the  ellipse  in  Grid  2 
had  a  =  -0.0032  and  a  =  0.0958  for  zero  angle  of 
attack.  This  variation  in  and  matches  Pulliam's 
findings  for  the  solution  of  the  Euler  equations  for 
subsonic  flow  past  ellipses  (21) ,  Pulliam  found  minor 
variations  in  the  grid  refinement  and  algorithm  type 
resulted  in  a  "general  lack  of  consistency  in  the  results" 
(21:1).  A  numerical  error  is  being  introduced,  however,  the 
specifics  are  not  clearly  understood  (22) . 

Figure  16  shows  the  Cp  for  both  grids  along  the  upper 
surface.  The  flow  conditions  are  Mach  =  0.3  and  a  =  0°. 
on  the  lower  surface  is  graphically  identical  with  the  upper 


61 


surface;  therefore  plotting  it  would  be  redundant.  The 


variation  in  Cp  at  the  nose  and  tail  between  Grid  1  and 
Grid  2  accounts  for  the  difference  in  C^.  The  variation  in 
can't  be  graphically  determined  from  the  Cp  plot. 


X 


Figure  16:  Cp  vs  X  (Upper  Surface) 

Figure  17  and  18  show  the  pressure  contours  for  Grid  1 
and  2  at  the  same  flow  conditions  as  in  Figure  16.  The 
contour  lines  in  both  graphs  are  scaled  to  the  same  minimum 
and  maximum  dimensionless  pressure,  0.6951  and  0.7558, 
respectively.  These  are  the  actual  minimum  and  maximum 
pressures  calculated  for  Grid  1. 

By  close  inspection,  one  can  see  the  higher  pressure 
along  the  upper  surface,  thus  creating  the  negative  lift. 
The  actual  minimum  and  maximum  pressures  for  Grid  2  are 
0.6939  and  0.7693.  The  pressure  range  for  Grid  2  is  larger 


than  that  of  Grid  1.  Pulliam  observed  similar  pressure 
contour  variations  when  different  grids  were  used  (21:1). 


Figure  17 :  Pressure  Contours  for  Grid  1 


Figure  18:  Pressure  Contours  for  Grid  2 


63 


V.2  Trajectory  Analysis 


The  trajectories  for  Grid  1  and  2  were  compared  at 
different  time-steps  and  with  different  Jacobi an-update 
methodologies.  The  initial  lift,  drag,  and  moment  on  the 
ellipse  was  calculated  from  the  steady-state  solution.  The 
ellipse  was  then  allowed  to  freefall  through  an  arbitrary 
trajectory  as  calculated  by  the  3-DOF  algorithm.  The  3-DOF 
algorithm  calculates  the  new  position  relative  to  an 
aircraft  flying  straight  and  level  at  constant  velocity. 

The  At  used  in  the  3-DOF  algorithm  for  both  grids  was 
calculated  by  using  a  Courant  number  of  10  and  Eg.  102. 

At  =  (A^) 

max|A,jt| 


The  minimum  time-step  for  Grid  1  was  At  =  0.0005  and 
At  =  0.0003  for  Grid  2.  The  varying  geometries  between  Grid 
1  and  2  caused  the  values  of  At  to  be  different.  The  3-DOF 
algorithm  used  the  C^,  and  moment  coefficient,  Cjjj, 
calculated  after  each  time-step  to  determine  the  new 
location  of  the  ellipse.  The  trajectory  calculation  was 
limited  to  1000  iterations  due  to  computer  resources.  This 
moved  the  ellipse  approximately  2  diameters  from  its 
original  location.  This  would  be  inadequate  in  moving  a 
store  out  of  the  influence  of  an  aircraft. 
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Three  different  cases  were  investigated,  each  with 
different  initial,  flow-field  conditions  and  ellipse  mass- 
properties.  Table  1  shows  these  three  different  cases. 


TABLE  1:  Initial  Flow  Conditions  and  Mass-Properties 


CASE: 

Mach 

X2 

Xcoq 

Mass* 

Iz** 

1 

0.3 

O 

o 

o 

1.0 

0.000 

1696.5 

475.0 

2 

0.3 

O 

O 

o 

0.3 

0.295 

3390.2 

1165.3 

3 

0.3 

-5.0° 

0.3 

0.295 

3390.2 

1165.3 

All  lengths  are  dimensionless 
*  (kg)  , 

**  (kg-m^) 

First,  the  trajectory  of  the  ellipse  COG  for  the 
potential-flow  solution  and  Grids  1  and  2  are  compared  using 
the  Case  1  properties.  At  =  0.0005  was  also  used  in  the 
potential-flow  solution.  Figure  19  shows  the  COG  dropping 
straight  down  for  the  potential-flow  solution  since,  the 
ellipse  generates  no  lift  or  drag.  The  small  error  between 
the  potential  flow  and  implicit  algorithm  COG  trajectories 
is  cause  by  the  numerical  error  in  the  implicit  algorithm, 
which  calculates  drag  on  the  symmetric  ellipse.  Grid  1  and 
2  show  close  agreement  through  the  1000  iteration  COG 
trajectory.  However,  the  original  assumption  was  to 
compared  the  pitch  angle  of  the  ellipses  between  the  two 
algorithms  and  not  the  x-  and  y-location.  Figure  20  shows 
the  pitch  angle,  0,  versus  the  y-COG  location. 
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Figure  19;  COG  Location  vs  Solution  Type  (Case  1) 


Figure  20:  0  vs  y-COG  (Case  1) 


The  general  trend  of  0  for  Grid  1  and  2  is  the  same.  T 
difference  in  magnitude  can  be  accounted  for  due  to  the 


different  initial  and  values  between  Grid  1  and  2. 

The  trend  for  the  potential-flow  solution  shows  a  different 
initial  slope  than  the  implicit  algorithm  results.  This 
variation  in  9  could  be  a  result  of  the  initial  properties 
of  the  flow-field  and  the  ellipse.  The  initial  conditions 
in  Case  1  are  "zero"  conditions,  where  no  lift,  drag,  or 
moment  is  generated  on  the  ellipse  from  the  potential-f lov/ 
algorithm.  Any  small  error  could  introduce  large  variations 
in  the  final  results. 

The  p-residual  jumps  from  the  steady-state  value  (10“®) 
to  a  value  on  the  order  of  10“^  due  to  the  initial  movement 
for  both  grids.  The  residual  then  gradually  decreases  down 
to  a  value  on  the  order  10”^  through  the  1000  iterations. 
This  trend  agrees  with  results  obtained  by  Simpson  (23) .  By 
observing  the  behavior  of  the  residual,  one  can  access  the 
likelihood  that  the  implicit  algorithm  is  generating  a 
reasonable  result.  As  long  as  the  residual  remains  small,  a 
value  on  the  order  of  10”^,  the  algorithm  results  should  be 
reasonable. 

Case  2  is  the  first  attempt  at  alleviating  the 
discrepancy  in  the  0-trend  between  the  potential  and 
implicit  algorithm  results.  Case  2  simply  moves  the  COG 
forward  of  the  center-of-pressure,  therefore  producing  a 
more  stable  store.  Figure  21  shows  the  0  versus  y-COG 
trajectory  for  the  Case  2  conditions. 

.Case  2  initial  0-trends  betv/een  the  potential  flow  and 
Grid  1  vary  the  same  as  they  do  in  Case  1.  The  initial 
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slope  of  the  0  versus  y-COG  for  the  potential-flow  result  is 
approximately  -0.001;  whereas  the  initial  slope  for  Grid  1 
is  approximately  -146.7.  The  "zero"  initial  conditions 
again  could  cause  the  variation  in  the  trends  between 
potential  flow  and  Grid  1. 


Case  3  is  an  attempt  at  removing  the  "zero"  initial 
conditions  from  the  problem.  In  this  case,  we  use  the  same 
store  mass-properties  as  in  Case  2,  but  start  the  ellipse  at 
an  a  =  -5°,  Figure  22  shows  the  variation  in  0  versus  y-COG 
through  a  1000  iterations.  The  ellipse  continues  to  pitch 
down  in  the  potential-flow  case,  while  the  ellipse  in  Grid  1 
is  trying  to  stabilize.  Case  3  indicates  there  are  large 
differences  between  the  potential  flow  and  Euler  implicit 
algorithm  results,  which  make  applicable  comparisons  very 
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Figure  22:  0  vs  Y-COG  (Case  3) 

V. 3  Implicit  Algorithm  Investigation 

The  effects  of  multiplying  the  minimum  time-step.  At, 
by  10  and  100  were  investigated  for  the  implicit  algorithm. 
The  purpose  is  to  evaluate  the  applicability  of  using  larger 
time-steps  in  the  implicit  algorithm  in  generating  a 
solution.  Figure  23  shows  the  COG  trajectory  results  and 
Figure  24  shows  0  versus  y-COG  for  Grid  1  with  Case  1 
initial  conditions. 

The  minimum  time-step.  At,  and  the  AtxlO  COG 
trajectories  compared  well  for  Grid  1,  however  the  AtxlOO 
results  are  extremely  poor.  The  minimum  time-step  may  be 
too  large  for  the  3-DOF  routine,  since  not  even  the  first 
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iteration  has  the  same  slope  as  the  other  results. 

Typically,  6-DOF  routines  (Reference  17)  utilize  a 
time-step  between  At  =  0,005  to  0.01  for  store  separation 
(26) .  The  time-step  for  the  AtxlOO  case  is  larger  than  a 
typical  6-DOF  time-step.  The  minimum  time-steps  of  0.0003 
and  0.0005  do  not  negatively  impact  the  3-DOF  routine, 
except  in  requiring  more  iterations  to  reach  the  solution. 

The  trend  for  0  in  Figure  24  matches  well  for  the  At 
and  AtxlO  results,  however  the  AtxlOO  results  are  extremely 
poor.  The  larger  period  and  amplitude  of  the  AtxlO  result 
compared  to  the  At  results,  can  be  accounted  due  to  the 
larger  time-step.  This  larger  time-step  is  than  multiplied 
by  the  lift  and  drag  values  calculated  by  the  implicit 
algorithm  solution  causing  the  larger  pitching  amplitude. 


Figure  23:  Grid  1  Trajectories  at  At,  AtxlO,  AtxlOO 
(Case  1) 
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Figure  24:  0  vs  y-COG  at  At,  AtxlO,  AtxlOO  for 

Grid  1  (Case  1) 

The  At  was  varied  further  to  investigate  its  affect  on 
the  trajectory.  Figure  25  shows  the  trajectories  for  Atx2, 
Atx4,  Atx6,  and  Atx8.  The  results  shows  that  the  magnitude 
of  the  At  directly  affects  the  magnitude  of  the  oscillation 
of  the  ellipse.  Since,  the  "correct”  trajectory  is  unknown, 
determining  a  valid  time-step  is  not  possible.  However,  the 
potential-flow  trajectory  of  AtxlO  matches  the  original  At 
trajectory,  as  shown  in  Figure  26.  This  is  an  indication 
that  the  3-DOF  algorithm  is  calculating  valid  solutions  at 
the  larger  time-steps. 

The  p-residual  trend  for  the  Atx2,  4,  6,  8,  and  10 
cases  matched  the  minimum  time-step  trend  addressed  earlier 
in  this  section.  The  p-residual  for  the  AtxlOO  case  jumped 
from  the  steady-state  case  to  a  value  on  the  order  10”^  and 
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remained  at  this  level.  This  is  an  indication  that  the 


AtxlOO  time-step  is  too  large  for  the  implicit  algorithm 
(23)  . 


Figure  25:  0  vs  y-COG  for  Varying  At's 

(Case  1) 


Figure  26:  6  vs  y-COG  for  the  Potential  Flow 

Solution  for  At  and  AtxlO  (Case  1) 
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The  results  of  the  Grid  1,  AtxlOO  case  indicated  the 
investigation  of  the  AtxlOO  case  for  Grid  2  would  not  be 
necessary  and  was  not  examined.  Figure  27  displays  COG 
trajectories  for  Grid  2,  Case  1  results;  there  is  good 
agreement  between  the  At  and  the  AtxlO  solutions.  Figure  28 
show  0  vs  y-COG  at  At  and  AtxlO  for  Grid  2,  Case  1,  The 
initial  0-trend  is  the  same  as  the  Grid  1  results.  The 
AtxlO  solution  allows  for  a  faster  determination  of  the  COG 
trajectory. 


Figure  27:  Grid  2  Trajectory  with  At  and  AtxlO 
(Case  1) 
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Figure  28:  9  vs  y-COG  at  At  and  AtxlO  for  Grid  2 

(Case  1) 

The  final  phase  of  the  investigation  involved  the 
examination  of  alternative  f lux-Jacobian  update 
methodologies.  In  work  done  by  Belk  and  Arabshabi 
(References  2  and  8) ,  the  flux-Jacobians  were  updated  after 
every  iteration  during  the  dynamic-grid  calculations. 
Simpson  found  it  wasn't  necessary  to  update  the  flux- 
Jacobians  after  every  time-step  (23) .  He  found  that  by 
calculating  the  flux-Jacobians  once,  they  could  be  used 
throughout  a  trajectory.  Simpson  had  also  unsuccessfully 
investigated  updating  the  flux-Jacobians  at  varying 
intervals  during  the  trajectory.  It  was  recommended  for 
this  research  the  flux-Jacobians  be  calculated  once  before 
the  trajectory  begins  and  frozen  through  the  trajectory. 

The  computational  time  decreased  by  approximately  45%  by 
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using  this  update  methodology.  Grid  1  and  2  were  evaluated 
with  the  flux-Jacobians  calculated  once  and  frozen 
throughout  a  trajectory  (1000  iterations) .  Figure  29 
compares  the  COG  trajectory  for  the  two,  Jacobian-update 
methodologies.  The  results  match  the  trends  found  by 
Simpson  (23) . 

As  an  additional  investigation,  the  Jacobian-update 
methodology  was  also  examined  for  the  AtxlO  time-step  for 
Grid  1  and  2.  Figure  30  shows  the  results  of  this 
comparison.  Figure  31  and  32  show  0  versus  y-COG;  the 
results  match  for  the  different  Jacobian-update 
methodologies.  The  increased  time-step  doesn't  affect  the 
trajectories  using  the  different  Jacobian  methodologies. 
These  results  suggest  the  alternative  Jacobian-update 
methodology  is  a  promising  one  and  should  be  investigated  in 
future  trajectory  analyses. 
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Figure  29:  COG  Trajectories  Utilizing  the 
Different  Jacobian  Update 
Methodologies  (Case  1) 


Figure  30:  COG  Trajectory  Utilizing  the 
Different  Jacobian-Update 
Methodologies,  AtxlO  (Case  1) 
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0  (RADIANS) 


Figure  31:  9  vs  y-COG  Utilizing  the  Different 

Jacobian  Update  Methodologies 
(Case  1) 


0  (RADIANS) 


Figure  32:  0  vs  y-COG  Utilizing  the  Different 

Jacobian  Update  Methodologies, 
AtxlO  (Case  1) 


Chapter  VI;  Conclusions  and  Recommendations 


VI. 1  Conclusions 

The  Grid  1  and  2  trajectories  matched  well  for  the 
ellipse;  the  0-trend  for  Grid  1  and  2  were  also  similar. 

This  demonstrates  the  potential  for  using  a  dynamic  grid, 
such  as  Grid  2  in  future  research.  The  numerical  error 
introduced  by  the  implicit  algorithm  made  comparisons  of  the 
Cog  trajectories  between  the  potential-flow  solution  and 
Grid  1  and  2  difficult. 

The  Case  1  and  2  examinations,  between  the  potential 
flow  and  implicit  algorithm  solutions  indicated  the  9-trend 
is  highly  sensitive  to  the  "zero"  initial  conditions. 

6-trend  results  for  Case  1  and  2  for  the  potential  flow  and 
implicit  algorithm  didn't  match  well.  Case  3  was  an  attempt 
at  examining  the  0-trend  at  non-zero,  initial  conditions. 

The  Case  3  results  showed  the  potential  flow  and  Grid  1 
results  still  didn't  correspond  for  the  same  initial 
conditions.  In  Case  3,  both  ellipses  started  to  pitch  nose- 
down,  but  the  Grid  1  ellipse  began  to  stabilize  through  the 
trajectory.  The  potential-flow  solution  drove  the  ellipse 
to  a  0  =  22.9°  with  no  indication  that  it  was  going  to 
stabilize.  The  potential-flow  solution  at  non-zero  initial 
conditions  requires  further  investigation. 

The  AtxlO  trajectory  results  for  Grids  1  and  2  shov;ed 
comparable  trends.  The  variation  in  pitch  amplitude  and 
frequency  is  caused  by  the  larger  time-step.  The  AtxlO 
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trajectories  need  to  be  examined  for  the  ellipse  used  in 
Case  2  and  3.  The  absence  of  experimental  data  in  this 
research  doesn't  allow  me  to  quantify  one  trajectory  better 
than  another. 

The  Jacobian-update  methodology  resulted  in  a  45% 
reduction  in  computational  time  for  the  calculation  of  a 
trajectory.  The  close  agreement  in  the  resulting  COG- 
trajectories  and  0  suggests  the  Jacobian-update  schedule 
requires  continued  examination  in  future  research. 

The  final  finding  of  this  work  is  the  criticality  of 
ensuring  the  time-step  utilized  falls  within  the  limits  of 
3-DOF  and  implicit  algorithms.  If  the  At  is  larger  than  the 
algorithms  can  manage  the  results  will  be  in  error. 

VI . 2  Recommendations 

Further  investigations  should  examine  the  use  of 
higher-order,  3-DOF  routines  in  determining  the  stores 
trajectory.  A  higher-order,  3-DOF  solution  may  provide 
closer  results  between  the  Grid  1  and  2  trajectories 
examined.  Other  3-DOF  algorithms  may  generate  closer 
results  between  the  potential-flow  solution  and  the  implicit 
algorithm. 

Perturbation  theory  could  be  used  to  account  for  the 
Mach  number  effects  in  the  potential-flow  solution.  These 
results  could  then  be  compared  with  results  from  the 
implicit  algorithm.  This  analysis  may  uncover  the  origin  of 
the  numerical  error  present  in  the  Euler  algorithm.  By 
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using  perturbation  theory  in  the  potential-flow  solution, 
the  resulting  trajectories  and  trends  may  compare  better 
with  implicit  algorithm  results. 

Additional  research  should  also  focus  on  comparing  the 
implicit  algorithm/ 6-DOF  routine  with  actual  wind  tunnel/6- 
DOF  trajectories.  The  lack  of  experimental  data  for  ellip.«e 
trajectories  and  the  inherit  problems  with  the  implicit- 
algorithm  results  for  the  ellipse  made  for  poor  comparisons. 
Wind-tunnel  data  and  the  6-DOF  routine  detailed  in  Reference 
17  would  be  an  excellent  comparison  tool  with  the  implicit 
algorithm.  Wind-tunnel  data  includes  Cj^,  C^,  and  data 
for  the  store  throughout  the  entire  flow-field.  These  data 
are  utilized  in  the  6-DOF  routine  and  could  be  utilized  in 
evaluating  algorithm  results. 

Other  grid-generation  programs  could  also  be 
investigated.  The  discontinuity  in  the  cells  along  the 
branch-cut  in  the  grid  could  be  eliminated  by  other  grid- 
generators.  By  not  starting  with  the  same  grid,  comparisons 
and  conclusions  on  the  results  are  difficult  to  make. 

Further  investigation  of  store  trajectories  isn't 
recommended  with  an  ellipse.  The  inherit  problems  of  the 
ellipse  highlighted  by  Pullium  (Reference  21)  and  in  this 
research,  suggests  it  is  a  difficult  tool  to  use  in 
trajectory  analysis.  Direct  comparisons  with  experimental 
data  is  required  to  give  valid  results  and  conclusions  on 
the  accuracy  and  applicability  of  the  CFD  algorithm. 
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APPENDIX  A:  Trapezoidal  Integration 


The  calculation  of  the  lift,  drag,  and  moment  on  the 
ellipse  is  performed  using  the  trapezoidal-integration 
technique.  This  method  is  used  in  the  potential-flow 
algorithm.  The  potential-flow  algorithm  calculates  the 
pressure  at  each  grid  point  along  the  surface  of  the 
ellipse.  The  surface  pressure  is  integrated  to  determine 
the  lift,  drag,  and  moment.  Trapezoidal  integration 
averages  the  surface  pressure  calculated  at  each  grid  point 
to  determine  the  pressure  over  the  interval,  using 

p.  ,  =  Jli - ^  (103) 

Jti  2 


The  X-  and  y-components  of  the  surface  pressure  are 
summed  in  the  body  axis  defined  in  Figure  6.  The  x-  and  y- 
components  of  force,  and  are  determined  using 


(104a) 


2=1 


ii,  =  E 


(104b) 


i=l 


The  lift  and  drag  on  the  ellipse  are  then  calculated 

using 


L  =  YjgCOsa  -  X^sina 


(105) 
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(106) 


D  =  Y^stna  +  Xj^cosa 

a  =  0  for  the  lift  and  drag  calculation  used  in  the  implicit 
algorithm. 

The  moment  calculation  is  slightly  different,  since  the 
COG  must  be  taken  into  account.  The  moment  is  positive  in 
the  counter-clockwise  direction.  It  is  calculated  using 

n-l 

^  ~  ~  ^ielx  ^  ^icLx^ycerii  ~  Xltl^)  (107) 

1=1 


where 


_  ^1.1  +  -^1 


_  Yi*!  *  Yi 


(108a) 


(108b) 


These  force  calculations  were  validated  using  the  zero 
lift  and  drag  requirements  for  potential  flow  with  no 
circulation.  The  moment  calculation  was  validated  using 
Eq.  7.  Table  2  shows  differences  between  the  calculate  and 
analytical  solutions  for  different  a.  The  values  shown  in 
Table  2  are  for  a  101x101  grid.  By  refining  the  grid,  the 
A's  decreased  to  values  on  the  order  of  10“^  for  a  1001x1001 
grid. 
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Table  2:  Delta  Force  vs  a 


a  (deg) 

AL 

AD 

AM 

O 

• 

o 

o 

1.56x10“^^ 

3.27x10“^^ 

7.59x10"^'^ 

10.0° 

5.72x10“’^ 

8.39x10“'^ 

8.15x10“^ 

O 

• 

o 

o 

1.71x10“”^ 

2.17x10“*^ 

6.83x10“^ 

OJ 

O 

• 

o 

o 

3. 10x10“^ 

7.58x10“'^ 

2.62x10"^ 

o 

• 

o 

o 

o.esxio""^ 

9.62x10"'^ 

1.24x10“^ 

The  results  in  Table  2  show  the  trapezoidal-integration 
technique  to  be  reasonably  accurate  compared  to  the 
analytical  results. 


APPENDIX  B;  Determination  of  Mach  Number 


The  Whitfield  algorithm  was  designed  to  work  at 
subsonic,  transonic,  and  supersonic  velocities.  This 
research  focuses  its  attention  on  the  subsonic  region. 
Evaluation  of  the  Cp  on  the  surface  of  the  ellipse 
determined  the  actual  Mach  number  used  to  perform  the  store 
separation  analyses.  These  values  were  also  compared  to  the 
potential-flow  solution.  A  201x201  grid  about  an  ellipse 
with  a  semi-major  axis  of  1.0  and  a  semi-minor  axis  of  0.2 
was  used  to  evaluate  the  Cp.  The  outer  radius  of  the  grid 
was  19.8  or  approximately  10  chords  from  the  center. 

The  algorithm  computes  Cp  on  the  surface  of  the  ellipse 
at  seven  different  Mach  numbers.  Simpson  recommended 
investigating  values  around  Mach  =  0.3,  since  the  algorithm 
has  been  successfully  tested  at  these  values  (23) .  Each 
solution  converged  until  a  residual  value  for  density  (p)  of 
order  10~®  was  obtained.  Table  3  shows  the  results  of  the 
seven  Mach  numbers  and  their  respective  values  of  Cj,,  C^^, 
and  Cjjj.  Figure  33  and  34  show  the  Cp  along  the  upper 
surface  of  the  ellipse  at  an  angle  of  attack  of  0.0.  The 
pressure  coefficient  curves  are  graphically  symmetrical  for 
the  upper  and  lower  surfaces. 
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Table  3:  Force  Coefficient  at  Various  Mach  Numbers 


Mach  Number: 

Cl 

Cd 

Cm 

0.10 

-0.0001 

0.1719 

0.0000 

0.15 

0.0000 

0.1188 

0.0000 

0.20 

0.0000 

0.0921 

0.0000 

0.25 

0.0000 

0.0762 

0.0000 

0.30 

0.0000 

0.0658 

0.0000 

0.35 

0.0000 

0.0586 

0.0000 

0.40 

0.0000 

0.0534 

0.0000 

Figure  33:  Cp  vs  X  (Upper  Surface) 
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Figure  34:  Cp  vs  X  (Upper  Surface) 


The  Cp  curve  at  M  =  0.40  is  the  closest  to  the 
potential-flow  values,  however  M  =  0.40  is  too  high  for  this 
research.  The  rest  of  the  research  used  a  Mach  number  of 
0.30.  The  Cp  curve  at  Mach  =  0.30  is  relatively  smooth  and 
compares  well  with  the  potential-flow  result.  More 
importantly,  this  Mach  number  is  low  enough  so 
compressibility  effects  can  still  be  considered  negligible. 
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.  scheme  is  employed  to  find  the  split-fluxes  and  the  Steger/Warming  flux-vector  method  is  used  to 
'  calculate  the  flux-Jacobians.  The  potential  flow  and  implicit  algorithm  are  combined  with  a  three- 
;  degree-of-freedom  algorithm  to  evaluate  the  planar,  freefell  trajectories  of  a  simple  store  shape.  The 
*  research  uses  two  different  grid-modification  techniques  in  the  implicit  algorithm  evaluation. 

;  Data  collected  for  both  grids  utilized  the  minimum  time-step  in  the  three-degree-of-freedom 

i  algorithm  for  a  Couranl  number  of  10.  Two  test  cases  involved  updating  the  flux-Jacobians  after 
I  every  time-step  and  only  once  during  every  1000  iterations.  The  effect  of  multiplying  the  minimum 
,  time-step  by  factors  of  2,  4,  6,  8,  10,  and  100  were  also  examined.  The  potential  flow  and  implicit 
!  algorithm  trajectories  did  not  compare  very  closely.  The  various  At  and  Jacobian-update  results 
matched  rather  closely. 
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